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We present a versatile high-level programming-language implementation of nonlinear topology 
optimization. Our implementation is based on the commercial software package Femlab, and it 
allows a wide range of optimization objectives to be dealt with easily. We exemplify our method 
by studies of steady-state Navier-Stokes flow problems, thus extending the work by Borrvall and 
Petersson on topology optimization of fluids in Stokes flow [Int. J. Num. Meth. Fluids 2003; 
41:77-107]. We analyze the physical aspects of the solutions and how they are affected by different 
D . parameters of the optimization algorithm. A complete example of our implementation is included 

pl«l . as Femlab code in an appendix. 
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INTRODUCTION 



The material distribution method in topology optimization was originally developed for stiffness design of mechanical 
structures but has now been extended to a multitude of design problems in structural mechanics as well as to 
optics and acoustics 0, 0, IE 01 • Recently Borrvall and Petersson introduced the method for fluids in Stokes flow Q . 
However, it is desirable to extend the method to fluids described in a full Navier-Stokes flow; a direction pioneered 
• i-h ■ by the work of Sigmund and Gersborg-Hansen |7], |8|, |9( . 

In the present work we present such an extension by introducing a versatile high-level programming-language 
implementation of nonlinear topology optimization, based on the commercial software package Femlab. It has a 
O • wider range of applicability than the Navier-Stokes problems studied here, and moreover it allows a wide range of 
l— ~~ ', optimization objectives to be dealt with easily. 

Extending the topology optimization method to new physical domains generally involves some rethinking of the 
design problem and some "trial and error" to determine suitable design objectives. It also requires the numerical 
\^ . analysis and implementation of the problem, e.g., using the finite element method (FEM). This process is accelerated 
a lot by using a high-level FEM library or package that allows different physical models to be joined and eases the 
tasks of geometry setup, mesh generation, and postprocessing. The disadvantage is that high-level packages tend to 
have rather complex data structure, not easily accessible to the user. This can complicate the actual implementation 
of the problem because the sensitivity analysis is traditionally formulated in a low-level manner. 

In this work we have used the commercial finite-element package Femlab both for the solution of the flow problem 
and for the sensitivity analysis required by the optimization algorithm. We show how this sensitivity analysis can 
be performed in a simple way that is almost independent of the particular physical problem studied. This approach 
proves even more useful for multi- field extensions, where the flow problem is coupled to, e.g., heat conduction, 
convection-diffusion of solutes, and deformation of elastic channel walls in valves and flow rectifiers |icj . 

The paper is organized as follows: In Sec. ^ we introduce the topology optimization method for fluids in Navier- 
O Stokes flow, and discuss the objective of designing fluidic devices or channel networks for which the power dissipation 
is minimized. In Sec. IHII we express the Navier-Stokes equations in a generic divergence form that allows them to be 
solved with Femlab. This form encompasses a wide range of physical problems. We also work out the sensitivity 
analysis for a class of integral-type optimization objectives in such a way that the built-in symbolic differentiation 
tools of Femlab can be exploited. In Sec. IIVI we present our two numerical examples that illustrates different 
aspects and problems to consider: The first example deals with designing a structure that can guide the flow in 
the reverse direction of an applied pressure drop. The general outcome of the optimization is an i'-shaped channel, 
but the example illustrates how the detailed structure depends on the choice of the parameters of the algorithm. 
The second example deals with a four terminal device where the fluidic channel design that minimizes the power 
dissipation shows a Reynolds number dependence. As the Reynolds number is increased a transition occurs between 
two topologically different solutions, and we discuss how the position of the transition depends on the choice of initial 
conditions. Finally in the appendix we include a transcript of our Femlab code required for solving the second 
numerical example. The code amounts to 111 lines - excluding the optimization algorithm that can be obtained by 
contacting K. Svanberg [Til [l3| . 
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II. TOPOLOGY OPTIMIZATION FOR NAVIER— STOKES FLOW IN STEADY STATE 



Although our high-level programming-language implementation is generally applicable we have chosen to start on 
the concrete level by treating the basic equations for our main example: the full steady-state Navier-Stokes flow 
problem for incompressible fluids. 

We consider a given computational domain Q with appropriate boundary conditions for the flow given on the 
domain boundary <9f2. The goal of the optimization is to distribute a certain amount of solid material inside f2 such 
that the material layout defines a fluidic device or channel network that is optimal with respect to some objective, 
formulated as a function of the variables, e.g., minimization of the power dissipated inside the domain. 

The basic principle in the material distribution method for topology optimization is to replace the original discrete 
design problem with a continuous one where the material density is allowed to vary continuously between solid and 
void 0. Thus in our flow problem we assume the design domain to be filled with some idealized porous material of 
spatially varying permeability. Solid wall and open channels then correspond to the limits of very low and very high 
permeability, respectively. 

In the final design there should preferably be no regions at intermediate permeability since otherwise it cannot be 
interpreted as a solution to the original discrete problem. Alternatively it may be possible to fabricate the device 
from polymeric materials such as PDMS that naturally have a finite permeability to the fluid [l4j ]. 



A. Governing equations for flow in idealized porous media 

We assume that the fluid flowing in the idealized porous medium is subject to a friction force f which is proportional 
to the fluid velocity v, c.f. Darcy's law. Thus f = -av, where a(r) is the inverse of the local permeability of the 
medium at position r. These properties of the idealized porous medium may only be approximately valid for an actual 
medium. However, the assumptions are not in conflict with any fundamental physical law, and since the converged 
solutions contain only solid walls and open channels, the specific nature of the idealized porous medium is of no 
consequence. 

The flow problem is described in terms of the fluid velocity field v(r) and pressure p(j). The governing equations 
are the steady state Navier-Stokes equation and the incompressibility constraint 

p(v ■ V)v = V • <r - av, (I) 
V • v = 0, (2) 

where p is the mass density of the fluid. For an incompressible Newtonian fluid the components <7ij of the Cauchy 
stress tensor er are given by 



/ dvi dvj 



x y =- P ^+^— + — ), (3) 



where r\ is the dynamic viscosity. The formalism is valid in three dimensions, but for simplicity we shall consider only 
two-dimensional problems, i.e., we assume translational invariance in the third dimension and set r = (x\,X2) and 
v = (t>i(r),i;2(r)). The boundary conditions will typically be either Dirichlet type specifying the velocity field v on 
the boundary or Neumann type specifying the external forces n • er. 

It is convenient to introduce a design variable field 7(1") controlling the local permeability of the medium. We let_7 
vary between zero and unity, with 7 = corresponding to solid material and 7 = 1 to no material. Following Ref. [lj 
we then relate the local inverse permeability a(r) to the design field 7(1") by the convex interpolation 

q\l - 7] 

a(7) = "min + ("max - "mm) ; , (4) 

q + i 

where q is a real and positive parameter used to tune the shape of a("f). Ideally, impermeable solid walls would be 
obtained with a max = oo, but for numerical reasons we need to choose a finite value for a max . For the minimal value 
we choose a m i n = 0.0 

For a given material distribution 7(1") there are two dimensionlcss numbers characterizing the flow, namely the 
Reynolds number 



(5) 
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describing the ratio between inertia and viscous forces, and the Darcy number 

Da = (6) 

describing the ratio between viscous and porous friction forces. Here £ is a characteristic length scale of the system 
and v a characteristic velocity. 

Almost impermeable solid material is obtained for very low Darcy numbers, in practice Da < 10~ 5 . Further insight 
into the meaning of the Darcy number is gained by considering Poiseuille flow in a channel or slit of width £ between 
two infinite parallel plate s of porou s materi al. In this case the fluid velocity inside the porous walls decays on a length 
scale £ Da , where £ Da = \f~Da£ = y/r)/a ma x.- See also Sec. IIV ATI for details on how the flow depends on Da. 



B. Power dissipation 

In the pioneering work by Borrvall and Petersson [| the main focus was on minimizing the power dissipation in the 
fluid. The total power $ dissipated inside the fluidic system (per unit length in the third dimension) is given by [TEj ] 



dr. (7a) 



In steady-state this is equal to the sum of the work done on the system by the external forces and the kinetic energy 
convected into it, 



<&(v,p,7)= / V" WiCTijVj - mvA^pvf) As. (7b) 
■>■>< > ^ - J 



Here n is a unit outward normal vector such that n • cr is the external force acting on the system boundary and n • er ■ v 
is the work done on the system by this force. Moreover, in the common case where the geometry and boundary 
conditions are such that the no-slip condition v = applies on all external solid walls, while on the inlet and outlet 
boundaries v is parallel to n and (n • V) v = 0,[2(j Eq. (|7b|l reduces to 

$(v,p,7)= / -n • v (p + \pv 2 ) As. (7c) 
J an 

Borrvall and Petersson showed that for Stokes flow with Dirichlet boundary conditions everywhere on the boundary 
d£l, the problem of minimizing the total power dissipation inside the fluidic device subject to a volume constraint on 
the material distribution is mathematically well-posed. Moreover it was proven that in the case where a("f) is a linear 
function, the optimal material distribution is fully discrete-valued. 

When a(-f) is not linear but convex then the solid/void interfaces in the optimal solution are not discrete zero/unity 
transitions but slightly smeared out. Convexity implies that the (negative value of the) slope of a at 7 = is 
larger than at 7 = 1; therefore there will be a neighbourhood around the discrete interface where it pays to move 
material from the solid side to the void. Using the interpolation in Eq. we have a'(0) = (a m in — awx)^ 2 and 
a'(l) — (a m i n — Q!max)x4^' ^ ar S e values of q the interpolation is almost linear and we expect almost discrete 
interfaces, whereas for small q we expect smeared out interfaces in the optimized solution. 

Consider the case when Eq. lt7c|l applies. If the system is driven with a prescribed flow rate then minimizing the 
total power dissipation is clearly equivalent to minimizing the pressure drop across the system. Conversely, if the 
system is driven at a prescribed pressure drop, then the natural design objective will be to maximize the flow rate 
which is equivalent to maximizing the dissipated power, c.f. Eq. (JTcJ. In either case the objective can be described 
as minimizing the hydraulic resistance of the system. 

For problems with more complex design objectives, such as a minimax problem for the flow rate through several 
different outlets, there will typically be no analog in terms of total dissipated power. In such cases there is no guarantee 
for the existence of a unique optimal solution and one has to be extra careful when formulating the design problem. 



III. GENERALIZED FORMULATION OF THE OPTIMIZATION PROBLEM 



For a given material distribution we solve the Navier-Stokes flow problem using the commercial finite element 
software Femlab. It provides both a graphical front-end and a library of high-level scripting tools based on the 



4 



Matlab programming language, and it allows the user to solve a wide range of physical problems by simply typing 
in the strong form of the governing equations as text expressions. The equations must then comply with a generic 
divergence form that eases the conversion to weak form required for the finite element solution. However, that is 
not a severe constraint since this is the natural way of expressing most partial differential equations originating from 
conservation laws. 

Since we have chosen fluidics as our main example, we begin by expressing the incompressible Navier-Stokes flow 
problem in divergence form. Then we state the optimization problem with a general form of the design objective 
function and perform the discretization and sensitivity analysis based on this generalized formulation. We stress 
that although for clarity our examples are formulated in two dimensions only, the method is fully applicable for 3D 
systems. 



A. The flow problem in divergence form 



Oil 


, r 2 = 


C12 


, r 3 = 


"0" 



021 




C22 





We first introduce the velocity-pressure vector u = [v±, v 2 , p] and define for i = 1,2,3 the quantities I\ and F t as 

(8) 

(9) 

(10a) 

(10b) 
(10c) 



and 

Fi = p(v ■ V)ui + a{j)vi, F 2 = p(v • V)v 2 + a(-f)v 2 , F3 = V • v. 
Using this, Eqs. and J5J) can be written in divergence form as 

V ■ Ti = Fi in 0, Governing equations 



Ri 
-n • r, 





Gi 



EdRj 



3=1 



dm 



on dtt, 
on dft, 



Dirichlet b.c. 
Neumann b.c. 



where Ti and Fi are understood to be functions of the solution u, its gradient Vu, and of the design variable 7. 
The quantity i2j(u, 7) in Eq. i|10b|) describes Dirichlet type boundary conditions. For example, fluid no-slip boundary 
conditions are obtained by defining R\ = v\ and R 2 = v 2 on the external solid walls. The quantity Gj(u, 7) in 
Eq. (|10c|) describe Neumann type boundary conditions, and [ii denote the Lagrange multiplier necessary to enforce 
the constraint Ri — 0, e.g., the force with which the solid wall has to act upon the fluid to enforce the no-slip 
boundary condition. Of course, it is not possible to enforce both Dirichlet and Neumann boundary conditions for the 
same variable simultaneously. Only when the variable Ui is not fixed by any of the Dirichlet constraints Rj does the 
Neumann condition Gi come into play, as all dRj/dui vanish and the Lagrange multipliers fj,j are decoupled from 
Eq. I|10cfl . Inactive Dirichlet constraints can be obtained simply by specifying the zero- function Ri = 0, that also 
satisfies Eq. <|10b|) trivially. 



B. The objective function 



In general the design objective for the optimization is stated as the minimization of a certain objective function 
$(u, 7). We shall consider a generic integral-type objective function of the form 

$(11,7)= /A(u, 7 )dr+ f J B(u, 7 )ds. (11) 
J n Jan 

In particular, we can treat the design objective of minimizing the power dissipation inside the fluidic domain by 
taking, c.f. Eq. (JgJ 



a = + -p-Y +J2 a M v f ' mQ and B = ° 011 dn - ( 12 ) 



dvi dvj \ 2 



\dxj dxi 



Alternatively, the objective of maximizing the flow out through a particular boundary segment dfl is obtained by 
choosing 

a n ■ ^ 1 n f — n • v on 90 o , 

A = m Q, and B = < n m ?L (13) 

[ on oil\oil , v ' 
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and objectives related to N discrete points r/. can be treated using Dirac delta functions as 

N 

A = J2 A k(^,l)S(r-r k ) in ft and 5 = on dfl. (14) 

k=l 

Finally we stress that not all optimization objectives lend themselves to be expressed in the form of Eq. I|ll(l - an 
example of which is the problem of maximizing the lowest vibrational eigenfrequency in structural mechanics. 

C. Optimization problem 

The optimal design problem can now be stated as a continuous constrained nonlinear optimization problem: 



min <&(!!, 7) (15a) 

7 

subject to : / 7(r)dr — /3|f2| < 0, Volume constraint (15b) 
Jn 

: < 7(1") < 1, Design variable bounds (15c) 

: Eqs. I|ll)a|) to (|ll)cl) . Governing equations (15d) 



With the volume constraint we require that at least a fraction 1 — ft of the total volume |fi| should be filled with 
porous material. 

The very reason for replacing the original discrete design problem with a continuous one by assuming a porous and 
permeable material, is that it allows the use of efficient mathematical programming methods for smooth problems. 
We have chosen the popular method of moving asymptotes (MMA) [llL Il2| , which is designed for problems with a 
large number of degrees-of- freedom and thus well-suited for topology optimization |3( • It is a gradient-based algorithm 
requiring information about the derivative with respect to 7 of both the objective function $ and the constraints. 
Notice that for any 7 the governing equations allow us to solve for u; therefore in effect they define u[7] as an implicit 
function. The gradient of $ is then obtained using the chain rule 

d r 1 r , si f d^> du , 

— $(117 ,7) = 0-+ / 7r-ir dr - 16 

L J 6*7 J n du 97 

However, because 11(7] is implicit, it is impractical to evaluate the derivative du/d~/ directly. Instead, we use the 
adjoint method to eliminate it from Eq. (|16fl by computing a set of Lagrange multipliers for Eqs. l|lUa|) to (|l(Jc|) 
considered as constraints 0|. For details see Sec. IIIIDl 

The optimization process is iterative and the fcth iteration consists of three steps: 

(i) Given a guess 7^) for the optimal material distribution we first solve Eqs. l|lUa|) to IjlUcll for u^ as a finite 
element problem using Femlab. 

(ii) Next, the sensitivity analysis is performed where the gradient of the objective and constraints with respect to 7 
is evaluated. In order to eliminate du/dj from Eq. I|16fl we solve the adjoint problem of Eqs. (|10afl to IjlOcfl for 
the Lagrange multipliers uW, also using Femlab. 

(iii) Finally, we use MMA to obtain a new guess ry( k+1 "> for the optimal design based on the gradient information 
and the past iteration history. 

Of the three steps, (i) is the most expensive computationalwise since it involves the solution of a nonlinear partial 
differential equation. 

D. Discretization and sensitivity analysis 

The starting point of the finite element analysis is to approximate the solution component u$ on a set of finite 
element basis functions {Vi,n( r )}> 



n 



(17) 



6 



where u,- in are the expansion coefficients. Similiarly, the design variable field 7(r) is expressed as 

70) = ^7n^4,«(r). 



(18) 



For our incompressible Navier-Stokes problem we use the standard Taylor-Hood clement pair with quadratic velocity 
approximation and linear pressure. For the design variable we have chosen the linear Lagrange element. |2l| 
The problem Eqs. QlOajl to I jlOcjl is discretized by the Galerkin method and takes the form 



L,(U, 7) - E N Ji A i = and M *( U ' = °' 

3=1 



(19) 



where Uj, Aj, and 7 are column vectors holding the expansion coefficients for the solution Uj, n , the Lagrange multi- 
pliers Hi t ni an d the design variable field 7„, respectively. The column vector contains the projection of Eq. Hll)a|) 
onto ifi_ n which upon partial integration is given by 



L7' 77 — 



/ [}Pi,n 

rj)dr+ / ipi. n Gids 
Jn Jon 



The column vector Mi contains the pointwise enforcement of the Dirichlct constraint Eq. (|10b(l 

M i>n = J R i (u(r i , n )). 



(20) 



(21) 



Finally, the matrix Ny = — dMi/dUj describes the coupling to the Lagrange multipliers in Eq. (|l(Jc|) . The solution 
of the nonlinear system in Eq. (|19(l above corresponds to step (i) in fcth iteration. The sensitivity analysis in step (ii) 
requires us to compute 



d<$> 



<9$ dVi 



(22a) 



which is done using the standard adjoint method |l6j |. By construction we have for any 7 that Lj(U(7), 7) — 
X}j=i Nj i A J (7) = and rVL„(U(7), 7) = 0. Therefore also the derivative of those quantities with respect to 7 is zero, 
and adding any multiple, say Uj and Aj, of them to Eq. I|22a|) does not change the result 



d 

dry 



$(U( 7 ), 7 ) 



l — l 



E 



3=1 



c?7 



,9M, 
97 



<9$ 



3=1 



' »=i 



OA 

Of 



-A (22b) 



Here we see that the derivatives dXJi /dj and dAi / d-y of the implicit functions can be eliminated by choosing U, and 
Ai such that 



E ( K ii U i ~ N Ji A i 



3=1 



and 



E N « G 3 

3 = 1 



(23) 



where we introduced Ky = —dhi/dTJj. This problem is the adjoint of Eq. (|19|) and U and A are the corresponding 
Lagrange multipliers. 

In deriving Eq. I|22b|) we implicitly assumed that N.y is independent of 7, i.e., that the constraint i?i(u, 7) is a 
linear function. If this is not true then the gradient d&/d~f computed from Eq. l|22b(l is not exact, which may leed 
to poor performance of the optimization algorithm if the constraints are strongly nonlinear. In order to avoid such 
problems it is necessary to include the nonlinear parts of the constraint vector M into L and move the corresponding 
Lagrange multipliers from A into U. While this is beyond the scope of the divergence form discussed in Sec. IIII AI it 
is certainly possible to deal with such problems in Femlab. Also the sensitivity analysis above remains valid since it 
relies only on the basic form of Eq. i|19[l for the discretized problem. 
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E. Implementation aspects 



We end this section by discussing a few issues on the implementation of topology optimization using Femlab. 

Firstly there is the question of how to represent the design variable 7(1"). The governing equations as expressed 
by r.j and F t in Eq. UlUafl depend not only on the solution u but also on 7, and the implementation should allow 
for this dependence in an efficient way. Here our simple and straightforward approach is to include 7 as an extra 
dependent variable on equal footing with the velocity field and pressure, i.e., we append it to the velocity-pressure 
vector, redefining u as 



[vi,V2,p,i\. 



(24) 



This was already anticipated when we denoted the basis set for 7 by {</>4.„(r)}. By making 7 available as a field 
variable we can take full advantage of all the symbolic differentiation, matrix, and postprocessing tools for analysing 
and displaying the material distribution. Appending 7 to the list of dependent variables we are required to define a 
fourth governing equation. However, since we are never actually going to solve this equation, but rather update 7 
based on the MMA step, we simply define 



F 4 = 0, G 4 = 0, i? 4 = 0. 



(25) 



It is crucial then that the finite clement solver allows different parts of the problem to be solved in a decoupled 
manner, i.e., it must be possible to solve Eqs. l|l(Ja|l - (|10c|) for ui for i = 1,2,3 while keeping Uj, i.e., 7, fixed. 

In Femlab the nonlinear problem Eq. (|19|) is solved using damped Newton iterations [TtJ ■ Therefore the matrices 
Ky = — dlii/dXJj and N,-j = — d'M.i/dXJj appearing in the adjoint problem Eq. (|23[1 are computed automatically as 
part of the solution process and can be obtained directly as Matlab sparse matrices. They are given by 



K 



m_ 



dT t 

dVu.j 



dr 



Jan ouj 



and 



N 



dRj 



(26) 



(27) 



Regarding the right-hand side vector d$/dXJi in Eq. J22J, notice that for a general objective as Eq. i|llfl . it has the 
form 



<9$ 



r dA 



dA 
dVu 



9-, 



dr 



on 



dB 

du. 



(28) 



It is not in the spirit of a high-level finite clement package to program the assembly of this vector by hand. In stead 
we employ the built-in assembly subroutine of Femlab. We construct a copy of the original problem sharing the 
geometry, finite element mesh, and degree-of-freedom numbering with the original. Only we replace the original fields 
IV Fi, and d with 



OA 
dVu- 



F 



dA 



and d 



dB 

du. t 



(29) 



Assembling the right-hand-side vector with this definition yields exactly Eq. Q28JI. c.f. Eq. J20J- An extra conve- 
nience in Femlab is that we can rely on the built-in symbolic differentiation tools to compute the derivatives dA/dui 
etc. In order to try out a new objective for the optimization problem, the user essentially only needs to change the 
text expressions defining the quantities A and B. 

After solving the adjoint problem Eq. (j23|l for U.; and to eliminate dXJi/d-y and dAi/dj for i = 1,2,3 in 
Eq. 122b|) we can evaluate the sensitivity 



d_ 



$(U, 7 ) 



d-y ^ V d 7 J J \ d-y 



1 A,- 



L 4 - J2 ( K J4 U ^ - N J 4 A, 



(30) 
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where = —dhi/d-f, N^4 = —dMi/dj, and L4 = d&/d~y in accordance with U4 = 7. Since the fourth variable 7 
is treated on equal footing with the other three variables, all expressions required to compute the matrices IQ 4 and 
N;.4 come out of the standard linearization of the problem. This is yet another advantage of including 7 as an extra 
dependent variable. 

When dealing with a problem with a volume constraint as in Eq. <|15b[) . it is necessary to compute the derivative 
of the constraint with respect to 7, 



d 



1 

W\ 



j(r) dr — (3 



1 

W\ 



(P„a(?) dr, 



(31) 



which can be obtained by assembling L4 with T4 = 0, F4, = 1, and G4 = 0. In the appendix we have included a 
transcript of the code required to set up and solve the example from Sec. IIV 51 below with Femlab. It amounts to 
111 lines of code, of which the majority are spent on setting up the actual Navier-Stokes flow problem. Only a minor 
part goes to set up the adjoint problem and perform the sensitivity analysis. Moreover, this part contains almost 
no reference to the actual physical problem being solved, and therefore it should apply for any multi-field problem 
expressed in the divergence form Eqs. (|10afl to (|10c|) with an objective function of the form of Eq. JUJ. The code 
example employs, but does not include, a Matlab implementation of the MMA optimization algorithm [Til IT^. IT^I] . 



1. Mesh dependence and regularization techniques 

It is well known that many topology optimization problems have trouble with mesh dependence. E.g. in stiffness 
design of mechanical structures it often pays to replace a thick beam with two thinner beams for a given amount of 
material. As the finite element mesh is refined, smaller and smaller features can be resolved and therefore appear 
in the optimized structure. In that sense the flow problem that we consider here is atypical because it is generally 
unfavorable to replace a wide channel with two narrower channels; hence the proof for the existence of a unique 
optimal solution with respect to minimization of the total power dissipation in Ref. 

The problem with mesh dependence can be overcome by various regularization techniques based on filtering of 
either the design variable 7(1*) or the sensitivity d^/d^f [3(- The regularization works by defining a certain length 
scale ro below which any features in 7(1") or d^/dj are smeared out by the filter; in both cases this results in optimized 
structures with a minimal feature size ~ ro independent of the mesh refinement. Unfortunately Femlab does not 
come with such a filter, and hence its implementation is an issue that has to be dealt with before our methodology 
here can be succesfully applied to problems that display mesh dependence. 

One strategy is to implement the convolution operation of the filter directly [3( . If the computational domain is 
rectangular and discretized by square finite elements this is both efficient and fairly easy to program, if not one simply 
uses a standard filter from the Matlab Image Processing Toolbox. For an unstructured mesh of triangular elements 
the programming is more involved and slow in Matlab due to the need to loop over the design variable nodes and 
searching the mesh for neighbouring nodes within the filter radius. Therefore an explicit matrix representation of the 
filter would often be preferred |isj . 

Another possible strategy is to solve an artificial diffusion problem for the design variable 7(1*) over some period 
in "time" At = r^/k where k is the "diffusion" constant. The diffusion equation could be included into the fields of 
Eq. I|29|) that are otherwise unused, and the "time" evolution solved using the built-in timestepper in Femlab. This 
procedure is equivalent to the action of a filter with Gaussian kernel of width ro , and it conserves the total amount of 
material during the filter action. The same approach could be used to smooth out the sensitivity. However, because 
d&/d"f n is sensitive to the local element size one would need to rescale it with J Q <^4 )Tl (r) dr before application of 
the filter - actually this is true for any filter acting on d^/dj whenever the mesh is irregular and J n </>4. n (r) dr not 
constant for all n. 

The major disadvantage of this strategy is that it involves solving a time evolution problem in each design iteration 
which could easily turn out to be the most time-consuming step. Alternatively the timestepping algorithm could be 
implemented by hand, e.g., deciding on the Crank- Nicholson algorithm with a fixed stepsize St < At. The mass and 
stiffness matrices for the diffusion problem can be obtained from Femlab, and the corresponding iteration matrix 
need only be factorized once for the given stepsize and could thus be reused in all subsequent design iterations, making 
this approach relatively cheap, although more cumbersome than using the built-in timestepper. 



2. Large-scale problems 



For large scale problems and three dimensional modeling it is often necessary to resort to iterative linear solvers 
because the memory requirements of a direct matrix factorization becomes prohibitive. In that case the strategy we 
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have outlined here of obtaining the K and N matrices directly as sparse matrices in Matlab and simply transposing 
K before the solution of the adjoint problem may not be practical. Alternatively, if the original physical problem is 
expressed in divergence form then the Femlab representation of that problem contains the symbolic derivatives of 
Ti, Fi, and Gi appearing in Eq. i|2ti|) . These fields can be transposed and set in the auxiliary copy of the original 
problem such that it effectively defines K,j = Kj^, while retaining the definitions in Eq. (|29|) for the right-hand-side 
vector Then the adjoint problem Eq. Il2.'ill can be solved without direct handling of the matrices in Matlab, and 
using the same iterative solver algorithm as would be employed for the original physical problem. Ultimately we still 
require an explicit representation of the matrices K^4 and N^4 to evaluate the sensitivity d^/d-y in Eq. (|30|) . 

From our point of view the major advantage of using Femlab in its present stage of development for topology 
optimization is not in solving large scale problems, though, but rather in the ease of implementation and the ability 
to handle problems with coupling between several physical processes. 

IV. NUMERICAL EXAMPLES 

In this section we present our results for topology optimization of Navier-Stokes flow for two particular model 
systems that we have studied. These systems have been chosen because they illustrate the dependence of the solution 
on the two dimensionless numbers Re and Da, measuring the importance of the inertia of the fluid and the permeability 
of the porous medium, respectively, relative to viscosity. Moreover we discuss the dependence of the solution on the 
initial condition for the material distribution. 

For simplicity and clarity we have chosen to consider only two-dimensional model systems. We note that the 
dimensionality of the problems has no fundamental consequence for the method and the numerics, but only affects 
computer memory requirements and the demand for CPU time. Our 2D examples can therefore be viewed as idealized 
test cases for our implementation of topology optimization. Yet, the 2D models are not entirely of academic interest 
only as they represent two limits of actual 3D systems. Due to planar process technology many contemporary lab- 
on-a-chip systems have a flat geometry with typical channel heights of about 10 /zm and widths of 1 mm, i.e., an 
aspect ratio of 1:100 0|. One limit is the case where the channel width is constant and the channel substrate and lid 
are patterned with a profile that is translation invariant in the transverse (width) direction. In the limit of infinitely 
wide channels the 2D-flow in the plane perpendicular to the width-direction is an exact solution, while it remains an 
excellent approximation in a 1:100 aspect ratio channel. This is the model system we have adopted for the numerical 
examples in the present work. The other important limit is when the channel width is not constant, but the channel 
height is sufficiently slowly varying that the vertical component of the fluid velocity can be neglected. Then writing 
the Navier-Stokes equation for the velocity averaged in the vertical (height) direction, the out-of-plane shear imposed 
by the channel substrate and lid gives rise to an absorption term -av. This approach was studied by Borrvall and 
Petersson [J, see also the footnote in Sec. Ill Al Thus, if one is interested in optimizing the height-averaged flow field 
in a flat channel the 2D model is sufficient. 

When solving the Navier-Stokes flow problem we use the standard direct linear solver in Femlab in the Newton 
iterations. Typically we have around 6000 elements in the mesh, corresponding to 30000 degrees-of-freedom. The 
constrained optimization problem is solved using a Matlab implementation of the MMA algorithm kindly provided 
by K. Svanberg except that we modified the code to use the globally convergent scheme described in Ref. [l2| . 

The example script included in the appendix employs only the basic algorithm mmasub, though. The design iterations 
are stopped when the maximal change in the design field is ||7( fc+1 ) — 7^- l || 00 < 0.01, at which point we typically have 

|$(fe+l) _ $(*)| < IO- 5 . 

A. Example: a channel with reverse flow 

Our first numerical example deals with the design of a structure that at a particular point inside a long straight 
channel can guide the flow in the opposite direction of the applied pressure drop. The corresponding problem with 
a prescribed flow rate was first suggested and investigated by A. Gersborg-Hansen We elaborate on it here to 
illustrate the importance of the choice of permeability for the porous medium. 

The computational domain is shown in Fig.^ It consists of a long straight channel of height £ and length L = 1Q£; 
the actual design domain, inside which the porous material is distributed, is limited to the central part of length b£. 
The boundary conditions prescribe a pressure drop of Ap from the inlet (left) to the outlet (right), and no-slip for 
the fluid on the channel side walls. 

The optimization problem is stated as a minimization of the horizontal fluid velocity at the point r* at the center 
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2.5£ 5£ 2.5£ 

FIG. 1: Computational domain for the reverse flow example. The design domain (gray) has length 5£ and height £, and the 
fluid enters and leaves the design domain through leads of length 2.5£. The boundary conditions prescribe a pressure drop of 
Ap across the system, and the design objective is to reverse the flow direction at the point r* at center of the channel. 



of the channel, i.e., the design objective is 

* = «i(r*). (32) 

In terms of the general objective Eq. (| 1 1 1) this is obtained with A = Ui(r)<5(r — r*) and B = 0. There is no explicit 
need for a volume constraint because neither of the extreme solutions of completely filled or empty can be optimal. 
When the design domain is completely filled with porous material we expect a flat flow profile with magnitude below 
Ap /(5£a max ). In the other extreme case when the channel is completely devoid of porous material the solution is 
simply a parabolic Poiseuille profile with maximum 

V A P 

V0 =8P—- (33) 

However, a structure that reverses the flow such that t>i(r*) becomes negative will be superior to both these extreme 
cases in the sense of minimizing $. 



1. Reverse flow m the Stokes limit, Re — 



We first consider the Stokes flow limit of small Ap where the inertial term becomes neglible. The problem is 
then linear and the solution is characterized by a single dimensionless parameter, namely the Darcy number Da, 
Eq. 10. We have solved the topology optimization problem for different values of Da. The initial condition for the 
material distribution was = 1, and the parameter q determining the shape of 01(7) in Eq. Q was set to q = 0.1. 
Anticipating that the structural details close to r* should be more important than those further away we chose a 
non-uniform finite element mesh with increased resolution around r* . 

Fig-Elshows the optimal structures obtained for Da = 10 -3 , 10~ 4 , 10~ 5 , and 10~ 6 . They all consist of two barriers 
defining an ^-shaped channel that guides the fluid in the reverse direction of the applied pressure drop. At Da = 10~ 3 
the two barriers are rather thick but leaky with with almost all the streamlines penetrating them; as the Darcy number 
is decreased the optimal structures become thinner and less penetrable. This result can be interpreted as a trade-off 
between having either thick barriers or wide channels. Thick barriers are necessary to force the fluid into the .S-turn, 
while at the same time the open channel should be as wide as possible in order to minimize the hydraulic resistance 
and maximize the fluid flow at the prescribed pressure drop. 

Notice that if we had chosen to prescribe the flow rate through the device rather than the pressure drop, then the 
optimal solution would have been somewhat different. When the flow rate is prescribed, it pays to make the gap 
between the barriers very small and the barriers very thick in order to force the fixed amount of fluid flow through 
the narrow contraction. The optimal structure is therefore one with a very large hydraulic resistance. In Ref. 8] this 
problem was circumvented by adding a constraint on the maximal power dissipation allowed at the given flow rate. 

In order to validate the optimality of the structures computed by the topology optimization we do as follows: For 
each of the optimized structures from Fig.|21we freeze the material distribution and solve the flow problem for a range 
of Darcy numbers. The resulting family of curves for ui(r*) vs. Da is shown in Fig. |3] where it is seen that each of 
the four structures from Fig. [21 do indeed perform better in minimizing «i(r*) than the others at the value of Da for 
which they are optimized. 

For Da < 10~ 5 the optimal value of v±(r*) tends to saturate because the thin barriers are then almost completely 
impermeable and the open channel cannot get much wider. In this limit the thickness of the optimized barrier 
structures approach the mesh resolution as seen in Fig.J^d). When the optimal barrier thickness gets below the mesh 
size we have observed the appearance of artificial local optima for the barrier structure. The problem is that the 
thin barriers cannot continuously deform into another position without going through an intermediate structure with 
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Da = 10" 4 




FIG. 2: Optimized structures (black) and streamlines at 5% intervals for Stokes flow (Re = 0) at Darcy numbers decreasing 
from 10 -3 to 10 -6 . Only the central part of length 2>i of the design domain is shown. The structures consist of two barriers 
defining an ^-shaped channel that reverses the flow at the central point r*. As the Darcy number is decreased, the optimized 
structures become thinner and less permeable. 
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FIG. 3: Comparing the performance of the structures from Fig.|5|optimized at fla op t for different values of Da. The objective 
«i(r*) is normalized with the velocity in an empty channel, vq, c.f. Eq. 1331 . 



barriers that are thicker by at least one mesh element. Depending on the initial condition, the optimization algorithm 
can therefore end up with a sub-optimal structure. We have tried to work around this problem by decreasing the 
value of q in order to smear out the solid/ void interfaces and thus reduce the cost of going through the intermediate 
structure. This did not work out well; the reason may be that the smearing property of a convex a("f) was derived 
for the objective of minimizing the power dissipation subject to a volume constraint. In the present example we are 
dealing with a different objective and have no volume constraint. However, when the barrier structures are resolved 
with at least a few elements across them the artificial local optima tend to be insignificant. Thus the problem can be 
avoided by choosing a sufficiently fine mesh, or by adaptively refining the mesh at the solid/void interfaces. 

Returning to Fig. [3]we notice that as Da increases all the structures perform poorly in minimizing v\ (r* ) , as they all 
approach vq . Extrapolating this trend one might suspect that the 5*-turn topology will cease to be optimal somewhere 
above Da = 1CP 3 simply because the porous material becomes too permable to make reversal of the flow direction 
possible. We have tested this hypothesis by performing an optimization at Da = 10 -2 , resulting in the structure 
shown in Fig. 0] where the value of the objective is fi(r*) = OTuo- It is seen to display a different topology from 
those of Fig. [21 with the design domain is almost completely filled with porous material blocking the flow through the 
channel. Only immediately above and below the point r* we see two empty regions emerging that act guide the flow 
away from r*. 

Actually, in all four cases from Fig. [21 starting from an empty channel the design iterations initially converge towards 
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Da = 10" 2 




FIG. 4: Optimized structure (black) and streamlines for Stokes flow at Da = 10 -2 ; only the central part of length 3£. The 
design domain is completely filled with porous material, except immediately above and below r* where two empty regions 
emerge. These voids divert the flow away from r*, resulting in a low velocity «i(r*) = O.lvo- 



(a) Ap = 0.2 x10 s , [fie = 23] (d) Ap = 0.2 x10 s , [fie = 26] 




(b) Ap = 0.5 x10 s , [fie = 42] (e) Ap = 0.5 x10 s , [fie = 47] 




(c) Ap = 1.0 x10 s , [fie = 64] (f) Ap = 1.0xl0 s , [fie = 71] 




FIG. 5: Optimized structures (black) and streamlines for Navier-Stokes flow; only a part of length 3.25^ near the center of 
the channel is shown. Panel (a)-(c) to the left show the optimized structures for different values of the control parameter 
Ap = Ap pi 2 J 'rf . For comparision the flow field when the optimized structure from Fig. EJc) is frozen and exposed to the 
elevated pressure drops is shown in panel (d)-(f) to the right. The Reynolds number is defined as Re = p£v max /r] where u max 
is the maximal velocity measured at the inlet; note that for a particular value of Ap, the Reynolds number is not fixed but 
differs slightly between left and right column. 

a symmetric structure blocking the flow like that in Fig.0] However at a certain point in the iterations an asymmetry 
in the horizontal plane is excited and the structure quickly changes to the two-barrier S'-geometry. Whether the 
optimization converge to an 5*- or an inverted S'-turn depends how the asymmetry is excited from numerical noise or 
irregularity in the finite element mesh; in fact the structure in Fig. E{b) originally came out as an inverted S but was 
mirrored by hand before plotting it to facilitate comparision with the three other structures. 

2. Reverse flow at finite Reynolds number 

We now consider flow at finite Reynolds number, characterized by the two dimensionless numbers Re and Da. The 
geometry and boundary conditions remain unchanged, for convenience we introduce a non-dimensional pressure drop 
Ap = App£ 2 /r] 2 , and finally we fix the Darcy number at Da — 10~ 5 . We note from Fig. [2] that this Darcy number 
allows some but not much fluid to penetrate the walls. We have nevertheless chosen this Darcy number for practical 
reasons, as the walls are "solid" enough and a lower value (more "solid" wall) would increase the calculation time. 

We have solved the topology optimization problem for different values of Ap, always using an empty channel as 
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21 L 21 

FIG. 6: Schematic illustration of the four-terminal device. Two inlet and two outlet leads (white areas) of height I and length 
21 are attached to the design domain (gray) of height 5£ and length L. The flow is characterized by the Reynolds number 
Re = p£v ma , x /ri, where u max is the maximal velocity at the inlets. 

initial condition. The results are shown in Fig.OJa)-(c) for Ap = 0.2, 0.5, and l.Ox 10 5 , where only a few streamlines are 
seen to penetrate the barriers. For comparision we also consider the flow field obtained when the structure optimized 
for Stokes flow at Da = 10~ 5 is frozen and exposed to the three different elevated pressure drops. This is shown in 
Fig. Eld)-(f): As the pressure drop is increased, more and more streamlines penetrate the barriers. Moreover we find 
a recirculation region emerging behind the second barrier which reduces the pressure drop over the neck between the 
barriers. 

Returning to Fig. |5fa)-(c), we find that the structures that have been optimized for the corresponding pressure 
drops are generally thicker than that optimized for Stokes flow, which reduces number of streamlines penetrating 
them. Also a beak-like tip grows on the second barrier that acts to bend the fluid stream down. Finally, on the back 
of the second barrier a wing- or spoiler-like structure appears that removes the recirculation. 

In summary, our first example has demonstrated that our implementation of topology optimization works, but that 
the optimal design and performance may depend strongly on the choice of the Darcy number. In particular, the zero 
Da limit solution contains zero thickness and yet impermeable barriers deflecting the fluid. In order to approximate 
this solution at finite Da and on a finite resolution mesh it is important to choose the Darcy number small enough 
that even thin barriers can be almost impermeable, but large enough to avoid difficulties with artificial local optima 
in the discretized problem when the barrier thickness decreases below the mesh resolution. 

B. Example: a four-terminal device 

Our second numerical example deals with minimization of the power dissipation in a four-terminal device subject 
to a volume constraint. The problem is found to exhibit a discrete change in optimal topology driven by the inertial 
term. The four-terminal device is related to one considered by Borrvall and Petersson for Stokes flow in Ref. Q ; the 
present example demonstrates that the optimization algorithm has difficulties in finding the optimal topology when 
there are two strong candidates for the global optimum. 

The computational domain, shown in Fig. consists of a rectangular design domain (gray) to which two inlet and 
two outlet leads (white) are attached symmetrically. The boundary conditions prescribe parabolic profiles for the flow 
at the inlets, zero pressure and normal flow at the outlets, and no-slip on all other external boundaries. Choosing the 
height I of the leads as our characteristic length scale, we define the Reynolds number as Re — plv max /r], where w max 
is the maximal velocity at the inlets. The Darcy number is fixed at Da = 10~ 4 to obtain reasonably small leakage 
through the porous walls. 

The optimization problem is stated as a minimization of the total power dissipation inside the computational 
domain, given by Eq. I|7a|l . subject to the constraint that at most a fraction (3 — 0.4 of the design domain should be 
without porous material, c.f. Eq. i|15b[) . 

Fig. E{a) and (b) shows the two optimal structures obtained for Re = 20 and 200, respectively, in a geometry with 
L = 3.5£. At Re = 20 the optimal structure turns out to be a pair of [/-turns connecting the inlets to the outlets 
on the same side of the design domain, while at Re = 200 the optimal structure is a pair of parallel channels. In 
order to minimize the power dissipation at low Re, the channel segments should be as short and as wide as possible, 
which favors the [/-turns in Fig.[7Ia). However, as the Reynolds number is increased, the cost of bending the fluid 
stream grows. When inertia dominates, larger velocity gradients appear in the long "outer lane" of the [/-turn. This 
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(a) Re = 20 (b) Re = 200 




FIG. 7: Optimal structures (black) and streamlines at 10% intervals for the four-terminal device at Reynolds number Re = 20 
and 200, respectively, in a geometry with L = 3.5^. 

increases the dissipation compared to low Re, where more fluid flows in the shorter "inner lane". At a certain point 
it will exceed the dissipation in the parallel channels solution 

*o = y (4+y) W 2 „ax, (34) 

as estimated from Poiseuille flow in two straight channels, each of length L + 4£ and height £. This number is 
independent of inertia due to translation symmetry, and we use $o as a natural unit of power dissipation (per unit 
length in the third dimension) in the following. 

Clearly the Reynolds number at which the transition between the two classes of solutions occurs will depend strongly 
on the ratio Ljl. For short lengths L <2£ the parallel channels solution is expected to be optimal at all Re, whereas 
for long lengths L > 3£ the ZJ-turn solution should be significantly better than the parallel channels solution at low 
Re. 



1. Dependence on the Reynolds number 

In the following we investigate more closely the transition between the [/-turns and the parallel channels solution 
as a function of the Reynolds number for the particular geometry L — it. The topology optimization problem is 
solved for different Re in the range to 200, using a homogeneous material distribution 7' ) = 0.4 as initial condition. 
For the parameter q determining the shape of 0(7) in Eq. Q we use a two-step solution procedure as suggested in 
Ref. 0| . First the problem is solved with q = 0.01 in order to obtain a solution with slightly smeared-out solid/void 
interfaces. Next this material distribution is used as initial guess for an optimization with q = 0.1 which generates 
fully discrete solid/void interfaces at the resolution of our finite element mesh. 

Fig.^a) shows the result for the normalized power dissipation $/$o obtained as a function of Re. At low Reynolds 
numbers the optimized solutions correctly come out as [/-turns with a power dissipation $ that is clearly less than 
$o- However, at high Reynolds numbers Re > 90 the method fails because the optimized solutions continue to come 
out as [/-turns even though this yields &/&o > 1- For Re > 160 the solution jumps from the simple [/-turns to a 
hybrid structure, as shown in the inset. The full lines in Fig. Ufa) show the result when the material distributions 
optimized for Re — 0, 50, and 180, respectively, are frozen and the power dissipation evaluated at different Re. It is 
seen that the optimized solutions, marked (o), all fall on or below the full lines which confirms that they are indeed 
superior to the other solutions of the [/-turn family. This also holds for Re > 90, except for the hybrid structures at 
Re > 160, that are actually inferior to the [/-turns. Moreover, at Re — 160 the optimized solution falls slightly above 
that optimized at Re = 180. This could be an indication that the hybrid structures are not local optima in design 
space after all, but rather a very narrow saddle point that the optimization algorithm has a hard time getting away 
from. 

The difficulty is that the two families of solutions, the [/-turns and the parallel channels, are both deep local minima 
for the power dissipation in design space. Using = 0.4 as initial condition, the initial permeability is everywhere 
very low, such that the porous friction almost completely dominates the inertia and viscous friction in the fluid. 
Therefore the iteration path in design space is biased towards low Reynolds numbers and the [/-turns solution. 

In order to circumvent this problem we have tried using a completely empty design domain with 7^ = 1 as initial 
condition. This should remove the bias towards the U -turns and allow the optimization algorithm to take inertia into 
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FIG. 8: Power dissipation $ in structures optimized for different Reynolds numbers; normalized with the Poiseuille flow 
result $o (dashed line), (a) Markers (o) show results when 7' ' = 0.4 is used as initial condition, failing to find the optimal 
solution for Re > 90. Full lines show the performance of the structures optimized at Re = 0, 50, and 180, when evaluated at 
different Reynolds numbers. As expected, all points fall on or below the full lines, except the hybrid solutions for Re > 160. 
(b) Comparision between the two different initial conditions 7^ = 0.4 (o) and 7' - 1 = 1 (x), showing the success of the empty 
channel initial condition in finding the optimal solution. The crosses (x) fall slightly below &/&o = 1 due to leakage through 
the porous walls (see the text). 
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FIG. 9: Comparision between structures optimized with the initially non-feasible material distribution 7' ' = 1 for different 
penalty parameters in the MMA optimization algorithm, revealing the difficulty in choosing the condition for finding the global 
optimum. Full line: the successful result from Fig. ISTbl with moderate penalty; (+): lower penalty yielding wrong result for 
40 < Re < 80; (□) higher penalty yielding wrong results for 80 < Re < 190. 



account from iteration one. The result is shown in Pig.jSfb). For Re < 80 the solutions are still ?7-turns, whereas for 
Re > 90 they come out as parallel channels. Notice that $/$o f° r the parallel channels solution is actually slightly 
less than unity, namely 0.98. This is due to a small amount of fluid seeping through the porous walls defining the 
device, which lowers the hydraulic resistance compared to the Poiseuille flow result derived for solid walls. 22] 

Strictly speaking the initial condition 7W) = 1 is not a feasible solution because it violates the volume constraint 
that at least a fraction 1 — j3 = 0.6 of the design domain should be filled with porous material. However, the MMA 
optimization algorithm penalizes this and reaches a feasible solution after a few iterations. This is controlled by 
choosing a penalty parameter. If the penalty for violating the constraint is small, the material is added slowly and 
only where it does not disturb the flow much. If the penalty is large, the material is added quickly and almost 
homogeneously until the constraint is satisfied. The succesful result from Fig. IHfb) was obtained with a moderate 
penalty. In Fig. |5| this is compared with results for smaller and larger penalty parameters, respectively. The figure 
shows that with the small penalization, the solution jumps to the parallel channels already at Re — 50 which is not 
optimal. For the large penalization, the solution does not jump until Re > 190. Also we observe hybrid structures 
similiar to those in Fig. |8{a) for Re > 130. We have thus not full control over the convergence towards the global 
optimum. 
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FIG. 10: Flow distributions at Re = 0, L = 3£ and q = 0.01. (a) Initial design field 7 (0) = 0.4. (b) Initial design field 
7 (0) = 0.4. (c) The optimal design field 7 W obtained at Da = 10~ 2 . 

2. Discussion of problems with local optima 

Further insight into the problem of local versus global optima is gained by inspecting the flow field in the initial 
material distribution 7 < -°- ) . This is shown in Fig. ^| for the Stokes flow limit, Re = 0. The streamlines are drawn 
as 10% contours of the streamfunction, and Fig. HOf a) shows that for <y(°> = 0.4 the streamline density is largest 
between the two leads on the same side of the design domain. Based on the sensitivity d^/dj the optimization 
algoritm therefore decides to remove material from these strong-flow regions in order to reduce the porous friction. 
The iteration path in design space is therefore biased towards the [/-turn solution. This remains true even at finite 
Reynolds numbers as long as the porous friction initially dominates inertia. 

Fig. HOf b) shows that when 7^°^ = 1 the streamline density is largest between the leads on the opposite side of the 
design domain. Because the volume constraint is violated the optimization algorithm has to place material somewhere, 
which it does in the weak-flow regions. The solution is therefore biased towards the parallel channels. Indeed if the 
penalty is chosen very small, the optimized solution comes out as parallel channels even for Stokes flow at L = 3£, 
which is far from optimal. When the penalty is larger and the material is added faster, we move away from this 
adiabatic solution and closer to the situation for 7W = 0.4. 

The additional complexity associated with making a proper choice of the penalty parameter is somewhat inconve- 
nient. We have therefore attempted to construct a more convex problem by increasing the initial permeability. This 
can be done either by increasing the Darcy number, or by decreasing the parameter q, c.f. Eq. Fig. HOf c) shows 
the optimal solution 7'*' obtained for Da = 10~ 2 and q = 0.01 at Re — 0. At this level the problem is convex because 
the solution is independent of the initial condition. Using this material distribution as initial guess and gradually 
decreasing the permeability to Da — 10~ 4 and 5 = 0.1 we correctly end up in the [/-turn solution. However, it is 
evident from Fig. HUf c) that 7**) has a fair amount of parallel channels nature. Using the same procedure of gradu- 
ally decreasing the permeability at higher Reynolds numbers therefore result in a transition to the parallel channels 
solution already for Re > 30, which is not optimal. Moreover, when the Reynolds number is increased and the inertia 
starts to play in, the system tends to loose convexity even at the initial high permeability. 

In summary, the topology optimization has difficulties in finding the global optimum for the problem. There are 
two strong candidates for the optimal structure, and the solution found is sensitive to the initial condition for the 
material distribution. Using an empty channel as initial condition, the method is able to find the correct solution for 
all Reynolds numbers. However, this successful result depends on a particular choice of the penalty parameter in the 
MMA algorithm. By using a high initial permeability of the porous medium, it is possible to convexify the problem 
at low Reynolds numbers, but continuation of this solution to the desired low permeability does not generally lead to 
the global minimum of the non-convex problem. 

In the original paper Ref. [lj it was argued that in Stokes flow the true optimal design should be rather insensitive 
to the choice of the Darcy number, although the dissipated power may deviate quite a lot from the zero Da limit. In 
our work we have observed that the actual solution found by the topology optimization may depend a great deal on 
the choice of the Darcy number, whereas the dissipated power should approach the zero Da limit roughly as \J Da. 

V. CONCLUSION 

Based on the work of Borrvall and Petersson we have extended the topology optimization of fluid networks to 
cover the full incompressible Navier-Stokes equations in steady-state. Our implementation of the method is based 
of the commercial finite element package Femlab, which reduces the programming effort required to a minimum. 
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Formulating the problem in terms of a general integral-type objective function and expressing the governing equations 
in divergence form makes the implementation very compact and transparent. Moreover the code for performing the 
sensitivity analysis should remain almost the same for any problem expressed in this way, whereas that required for 
describing the physical problem of course changes. Topology optimization of multi-field problems can therefore be 
dealt with almost as easy as a single realization of the underlying physical problem. 

We would like to mention that our methodology is not as such restricted to the (large) class of physical problems 
that can be expressed in divergence form. Femlab also allows problems to be stated directly in weak form, e.g., 
for systems with dynamics at the boundaries. This does in fact not invalidate the sensitivity analysis worked out in 
Sec. IIII Dl since this analysis only relies on the basic structure of the discretized nonlinear problem and the availability 
of the Jacobian matrix. It is therefore possible to apply our methodology to even larger classes of physical problems 
than the ones comprised by the divergence form. 

Our implementation of topology optimization has been tested on two fluidics examples in 2D, both illustrating the 
influence of different quantities and conditions on the efficiency of the optimization method. 

The first example, a channel with reversed flow, illustrates the influence of the Reynolds number Re and the Darcy 
number Da on the solutions. We have shown that the choice of Da has a strong impact on the solution when the 
structure contains barriers to deflect the fluid stream. 

The second example, minimization of the power dissipation in a four-terminal device, reveals the problems of 
determining the global minimum when two strong minima are competing. This problem is highly non-convex, and 
we have shown that the solution depends on the initial condition. For an initial homogeneous material distribution, 
the porous friction dominates and the solution does not come out as the global optimum in all cases. Using an empty 
channel as the initial state, inertia plays a role from the beginning, and better results can be obtained. However, 
this initial condition in fact violates the volume constraint, and the part of the optimization routine correcting this 
depends on a penalty factor. Unfortunately, the particular value chosen for this factor strongly influences the results. 
Increasing the Darcy number makes the problem more convex, but continuation from large to small Da, i.e., from 
high to low permeability of the porous material, does not generally end up in the global optimum. 

In conclusion, we have shown that our implementation of topology optimization is a useful tool for designing fluidic 
devices. 
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APPENDIX 

•/, FEMLAB CODE FOR THE 4-TERMINAL DEVICE EXAMPLE OF SEC. 4.2 
clear fern femadj 

DEFINE REYNOLDS NUMBER, DARCY NUMBER, LENGTH OF DESIGN DOMAIN, AND VOLUME FRACTION 
Re = 50; 
Da = le-4; 
L0 = 3.0; 
beta = 0.4; 

7. DEFINE GEOMETRY, MESH, AND SUBD0MAIN/B0UNDARY GROUPS [SEE FIG. 6] 

fem.geom = rect2(0,L0,0,5) + rect2(-2 ,0 , 1 ,2) + rect2(-2,0,3,4) + rect2(L0 ,L0+2, 1 , 2) ... 

+ rect2(L0,L0+2,3,4) ; 
fern. mesh = meshinit(fem, 'Hmaxsub' , [3 0.125]); 
7. subdomain groups 1: design domain 2 : inlet/outlet leads 
fem.equ.ind = {[3] [1 2 4 5]}; 

7. boundary groups 1: walls 2: inlets 3: outlets 4: interior 

fem.bnd.ind = {[2:3 5:8 10 12:14 16:18 20:22] [4 23] [1 24] [9 11 15 19]}; 

7. DEFINE SPACE COORDINATES, DEPENDENT VARIABLES, AND SHAPE FUNCTIONS 

fem.sdim = {'x' 'y'}; 

fern. dim = {'u' 'v' 'p' 'gamma'}; 

fern, shape = [2211]; 

7. DEFINE CONSTANTS 

fern. const . rho = 1; 

fern. const . eta = 1; 

fern. const .umax = Re; 

fern. const . alphamin = 0; 



fern. const . alphamax = 1/Da; 
fern. const . q = 0.1; 

PhlO = 96*f em. const. eta* (L0+4)*f em. const. umax~2/9; 

•/. DEFINE EXPRESSIONS ON SUBDOMAIN AND BOUNDARY GROUPS 

fem.equ.expr = {'A' 'eta*(2*ux*ux+2*vy*vy+(uy+vx)*(uy+vx))+alpha*(u*u+v*v) ' ... 

'alpha' { ' alphamin+(alphamax-alphamin) *q* (1-gamma) / (q+gamma) ' '0'}}; 
fem.bnd.expr = {'B' '0'}; 

•/. DEFINE GOVERNING EQUATIONS AND INITIAL CONDITIONS [SEE EQS. (8) AND (9)] 
fern. form = 'general'; 

fern. equ. shape = {[1:4] [1:3]}; '/, only define gamma on subdomaln group 1 

fem.equ.ga = {{{'-p+2*eta*ux' 'eta*(uy+vx) '} {'eta*(uy+vx) ' '-p+2*eta*vy'} {0 0} {0 0}}}; 
fern. equ. f = {{'rho*(u*ux+v*uy)+alpha*u' 'rho*(u*vx+v*vy)+alpha*v' 'ux+vy' 1}}; 
fern. equ. inlt = {{0 beta}}; 
•/. DEFINE BOUNDARY CONDITIONS 

f em. bnd. shape = {[1:3]}; '/, do not define gamma on any boundaries 

fem.bnd.r = {{'u' 'v' 0}... '/, walls: no-slip 

{'u*nx+4*umax*s*(l-s) ' 'v' 0}... '/, inlets: parabolic profile 

{0 'v' 0} ... '/, outlets: normal flow 

{0 0}}; '/, interior: nothing 

fern. bnd. g = {{0 0}}; '/, zero prescribed external forces everywhere 

•/. PERFORM LINEARIZATION, DEGREE-OF-FREEDOM ASSIGNMENT, AND ASSEMBLE INITIAL CONDITION 
fern = f emdif f (f em) ; 
fem.xmesh = meshextend(f em) ; 
fern. sol = asseminit (f em) ; 

•/. DEFINE STRUCTURE FOR COMPUTING RIGHT-HAND-SIDE IN ADJOINT PROBLEM [SEE EQ. (29)] 
femadj = fern; 

femadj .equ.ga = {{{ 'dif f (A,ux) ' 'dif f (A,uy) ' } { 'dif f (A, vx) ' 'dif f (A,vy) '} ... 

{'diff (A,px) ' 'diff (A,py) '} {'dif f (A.gammax) ' 'dif f (A.gammay) '}}}; 
femadj. equ. f = {{'diff (A,u) ' 'diff(A,v)' 'diff(A,p)> 'dif f (A, gamma) '}}; 
femadj. bnd. g = {{'dif f (B,u) ' 'diff(B,v)' >diff(B,p)> 'dif f (B, gamma) '}} ; 
femadj. xmesh = meshextend(f emadj ) ; 

•/. GET INDICES OF DESIGN VARIABLE IN THE GLOBAL SOLUTION VECTOR (fern. sol. u) 

i4 = f ind(asseminit (f em, 'Init' , {'gamma' 1} , ' Out ' , 'U' ) ) ; 

•/. COMPUTE VOLUME BELOW DESIGN VARIABLE BASIS FUNCTIONS 

L = assemble(f em, 'Out' ,{'L'}) ; 

Vgamma = L(i4) ; 

Vdomain = sum (Vgamma) ; 

7, GET INDICES OF VELOCITY-PRESSURE VARIABLES 

il23 = f ind(asseminit(f em, 'Init' ,{'u' 1 'v' 1 'p' 1} , 'Out ' , 'U' ) ) ; 

•/. DEFINE VARIABLES AND PARAMETERS FOR MMA OPTIMIZATION ALGORITHM [SEE REF. [11,12,13]] 

aO = 1; 

a = 0; 

c = 20; 

d = 0; 

xmin = 0; 

xmax = 1 ; 

xold = f em. sol .u(i4) ; 
xolder = xold; 
low = 0; 
upp = 1; 

7, DESIGN LOOP FOR THE ACTUAL TOPOLOGY OPTIMIZATION 
for iter = 1:100 

7. SOLVE NAVIER-STOKES FLOW PROBLEM TO UPDATE VELOCITY AND PRESSURE 
fem.sol = f emnlin(f em, 'Solcomp' ,{'u' 'v' 'p'} , 'U' , fern. sol .u) ; 
•/, SOLVE ADJOINT PROBLEM FOR LAGRANGE MULTIPLIERS 

[K N] = assemble(fem, 'Out' ,{'K' 'N'} , 'U' , fern. sol .u) ; 

[L M] = assembled emadj , 'Out' ,{'L' 'M'} , 'U' , fern, sol .u) ; 

femadj. sol = f emlin( ' In' , { 'K' K(il23, il23) ' 'L' L(il23) 'M' zeros (size (M) ) 'N' N(:,il23) 
7, SENSITIVITY ANALYSIS 
gamma = fern. sol .u(i4) ; 

Phi = postint(f em, 'A' , 'Edim' ,2) + postint(f em, 'B' , 'Edim' , 1) ; 
dPhidgamma = L(i4) - K(il23, i4) ' *f emadj . sol .u; 
7. PERFORM MMA STEP TO UPDATE DESIGN FIELD 
x = gamma; 

f = Phi/PhiO; g = gamma ' *Vgamma/Vdomain - beta; 

dfdx = dPhidgamma/PhiO ; dgdx = Vgamma ' /Vdomain ; 
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d2fdx2 = zeros (size (gamma) ) ; d2gdx2 = zeros (size (gamma' )) ; 

[xnew,y ,z, lambda, ksi , eta,mu,zeta, s , low,upp] = mmasub(l , length (gamma) , iter , ... 

x,xmin,xmax,xold,xolder , f ,df dx, d2f dx2 ,g,dgdx , d2gdx2 , low.upp , aO, a, c ,d) ; 
xolder = xold; xold = x; gamma = xnew; 
7. TEST CONVERGENCE 

if iter >= 100 I max(abs (gamma-xold) ) < 0.01 
break 

end 

7, UPDATE DESIGN VARIABLE 

uO = fem.sol.u; u0(i4) = gamma; 

fern. sol = femsol(uO); 

7, DISPLAY RESULTS FOR EACH ITERATION STEP 

disp(sprintf ('Iter. :°/,3d Obj . : '/.8.4f Vol.: °/,6.3f Change: 7,6. 3f, ... 

iter , f , xold' *Vgamma, max (abs(xnew-xold) ) ) ) 
postplot (f em, ' arrowdata' , { 'u' ' v' } , 'tridata' , 'gamma' , ' trimap' , 'gray ' ) 
axis equal; shg; pause (0.1) 

end 
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